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SUMMARY 

The current work is initiated in an effort to obtain an efficient, accurate, and robust algorithm 
for the numerical solution of the incompressible Navier-Stokes equations in two- and three- 
dimensional generalized curvilinear coordinates for both steady-state and time-dependent flow 
problems. This is accomplished with the use of the method of artificial compressibility and a 
high-order flux-difference splitting technique for the differencing of the convective terms. Time 
accuracy is obtained in the numerical solutions by subiterating the equations in pseudo-time 
for each physical time step. The system of equations is solved with a line-relaxation scheme 
which allows the use of very large pseudo-time steps leading to fast convergence for steady- 
state problems as well as for the subiterations of time-dependent problems. Numerous laminar 
test flow problems are computed and presented with a comparison against analytically known 
solutions or experimental results. These include the flow in a driven cavity, the flow over 
a backward-facing step, the steady and unsteady flow over a circular cylinder, flow over an 
oscillating plate, flow through a one-dimensional inviscid channel with oscillating back pressure, 
the steady-state flow through a square duct with a 90° bend, and the flow through an artificial 
heart configuration with moving boundaries. An adequate comparison with the analytical or 
experimental results is obtained in all cases. Numerical comparisons of the upwind differencing 
with central differencing plus artificial dissipation indicates that the upwind differencing provides 
a much more robust algorithm, which requires significantly less computing time. The time- 
dependent problems require on the order of 10 to 20 subiterations, indicating that the elliptical 
nature of the problem does require a substantial amount of computing effort. 
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PREFACE 


The use of computational fluid dynamics to analyze incompressible or very nearly incom- 
pressible flows has grown considerably in the last several years. The fact that the incompressible 
Navier-Stokes equations can be used as a governing approximation in a wide range of disciplines 
has spawned a need for a generalized flow solver. Generalized in the sense that the entire for- 
mulation for the basic code be independent of geometrical approximations. Besides being able 
to reliably produce trustworthy results, a useful analysis code should not require many different 
numerical parameters which have to be correctly specified before any answer can be obtained. 
As an analysis tool such a codes value can be measured by it 3 robustness and its ability to be 
used easily by people without extensive experience in running the code. 

The current work is an effort to improve upon an already successful incompressible flow 
solver known as INS3D. This code uses a very well known and successful approach, namely 
that of artificial compressibility, coupled with central differencing plus artificial dissipation, 
and solved with an approximate factorization scheme. This code, being limited to steady-state 
problems, also requires the specification of artificial dissipation coefficients, a time-step size, and 
an artificial compressibility parameter. The result of the current work is a code which takes the 
artificial compressibility approach a step further, by analyzing the physics of this approach, and 
applying a more appropriate differencing scheme. With the use of a line-relaxation solution, 
only one numerical parameter need now be specified, the artificial compressibility constant. 
The following work is intended to show that the resulting code is a robust means of obtaining 
steady-state and time-dependent solutions to the incompressible Navier-Stokes equations. 

This work was partially funded by NASA Marshall Space Flight Center. In addition, some 
of the funding for the artificial heart work came from the NASA Technology Utilization Office. 
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CHAPTER ONE 


INTRODUCTION 

1.1 Background and Motivation 

Numerical solutions to the incompressible Navier-Stokes equations are in greater 
demand than ever before as the field of computational fluid dynamics (CFD) increases 
its impact as an engineering tool. Problems which can be addressed by the incom- 
pressible Navier-Stokes equations include low-speed flows in aerodynamics, internal 
flows in propulsion, liquid flows in hydrodynamics, and problems in biomedical fluid 
analysis. With this wide range of disciplines come a wide range of possible problem 
geometries. Thus a useful flow solver will require that no simplifying assumptions be 
made about the geometries, and that the code easily accept different model configura- 
tions. This second criteria will require the use of generalized curvilinear (body-fitted) 
coordinates in the code formulation. The more efficient and robust that a code can 
become, the more useful a tool it will be for analysis. Therefore, there is a continuing 
interest in finding solution methodologies which will produce results using the least 
amount of computing time and with the least amount of effort by the investigator. 
This is particularly true for unsteady problems. Time-accurate solutions of the in- 
compressible Navier-Stokes equations are particularly time consuming because of the 
elliptical nature of the governing equations. A disturbance at one point in space affects 
the entire flow domain instantaneously. This requires that the numerical algorithm 
propagates information through the entire flow domain during one discrete time step. 
Thus some type of iterative scheme is usually required to solve the equations in time. 
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Another important item which directs the development of software is the current trend 
of new hardware availability. The newest generation of supercomputers has provided 
an order of magnitude increase in the available processor memory. This can be used 
to efficiently implement memory intensive algorithms which would otherwise be too 
costly. 

1.2 Available Solution Methods 


The primitive variable formulation of the incompressible Navier- Stokes equations 
is one in which the system of equations is written with the pressure and velocity 
components as the dependent variables. This form is perhaps the most widely used as 
a starting point for the development of a numerical solution method. Several methods 
do exist in which the incompressible Navier-Stokes equations are formulated into a 
nonprimitive- variable form. These include a velocity-vorticity method as proposed 
by Fasel [1], and a vector potential- vorticity method proposed by Aziz and Heliums 
[2]. Each of these methods require the solution of three Poisson equations at each 
time level when solving the three-dimensional (3-D) equations, which can become very 
computationally expensive. Further work on the velocity-vorticity method and the use 
of a direct solver is being implemented by Hafez et al [3]. These methods currently 
appear to be too expensive for solving large 3D problems. 

This leaves methods formulated in primitive variables. Most methods using prim- 
itive variables can be classified into three groups. The first of these, and historically 
one of the most commonly used primitive variable schemes, is the pressure Poisson 
method as first introduced by Harlow and Welch [4]. In this method, the velocity 
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field is advanced in time using the momentum equations. Then a Poisson equation 
in pressure, which is formed by taking the divergence of the momentum equations, is 
solved for the pressure at the current time level such that the continuity equation will 
be satisfied at the next time level. This method can be very costly because of the need 
to solve the Poisson equation at every time iteration, and because the velocity and 
pressure fields are only indirectly coupled. 

The second group of primitive- variable methods which has been used extensively 
is that of the fractional-step method originally introduced by Chorin [5], and used by 
Kim and Moin [6], and Rosenfeld et al. [7]. This method advances the solution in 
time using two (or more) steps. First the momentum equations are solved for an 
intermediate velocity field which will generally not be divergence free. In the next 
step, a pressure field is obtained which will map the intermediate velocity field into a 
divergence free field, thus obtaining the solution at the next time level. The second 
step generally requires the solution of a Poisson equation in pressure. This primitive 
variable method is similar to the first, and although it relates the changes in the 
pressure field more directly to the divergence of velocity, the pressure and velocity 
field solutions are still only indirectly coupled. 

The third primitive-variable method is that known as the artificial compress- 
ibility approach and was first introduced by Chorin [8]. It has been used with much 
success by Kwak et al. [9] for solving complex incompressible flow problems in gener- 
alized coordinates. In this formulation, a time derivative of pressure is added to the 
continuity equation. Together with the momentum equations, this forms a hyperbolic 
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system of equations which can be marched in pseudo-time to a steady-state solution. 
The method can also be extended to solve time-dependent problems [10] by using 
subiterations in pseudo-time at every physical time step. If all that is desired is the 
steady-state solution the artificial compressibility method can be a very efficient for- 
mulation because it does not require that a divergence free velocity field be obtained at 
each iteration. The artificial mechanism by which this corrects the flow field to satisfy 
the continuity equation is easy to understand: if in a computational cell the net flux 
of mass becomes greater than zero, thus making the divergence of velocity positive, 
the pressure in that cell is decreased and through the action of the pressure gradient 
increases the force drawing the fluid toward the cell; conversely if the divergence of 
velocity is negative, the pressure increases, which increases the pressure force pushing 
the fluid away. In this way the pressure and velocity fields are directly coupled. The 
addition of the time derivative of pressure to the continuity creates a hyperbolic sys- 
tem of equations complete with artificial pressure waves of finite speed. Since this is 
the case, many of the well-developed compressible flow algorithms can be utilized for 
this method. 


1.3 Objective and Approach 

The goal of the current work was to develop a method to solve the incompressible 
Navier-Stokes equations for both steady-state and time-dependent problems. This 
must be done as efficiently as possible so it can be used in solving complex 3-D problems 
in generalized coordinates which may require a large number (100,000 to 500,000) of 
grid points. In an effort to develop the methodology for such a code, work was done 
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in both two dimensions (2-D) and 3-D. The reason for developing both a 2-D and a 
3-D code is that many aspects of the code can be tested and developed more efficiently 
in 2-D, particularly when studying time-dependent problems. In addition, while most 
problems of interest are 3-D in nature, studying a 2-D model of the problem can often 
provide insight into the physics of the problem. 

The artificial compressibility approach was chosen for the current work because 
of the direct coupling between the velocity and pressure fields and because of the hy- 
perbolic nature of the resulting equations. One previous application of this method by 
Kwak et a/.[9] employed the use of central differencing of the convective fluxes, with 
the addition of artificial dissipation to prevent odd-even decoupling of the grid points. 
This method requires the specification of artificial dissipation coefficients which can 
strongly influence the resulting solution, and tends to degrade the accuracy of the 
solution [11]. In an effort to produce a code with as few numerical parameters as 
possible, upwind differencing was chosen to difference the convective fluxes. An ad- 
ditional advantage of upwind differencing when using an implicit scheme is that it 
will generally produce a more diagonally dominant system of equations than central 
differencing. Since the equations to be solved are hyperbolic, some of the upwind 
differencing schemes which have recently been developed for the compressible Euler 
and Navier-Stokes equations by a number of authors [12-15] can be utilized. Using 
the method of Roe [12] the convective terms are differenced by an upwind method 
that is biased by the signs of the eigenvalues of the local flux Jacobian. This is ac- 
complished by casting the governing equations in their characteristic form and then 
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forming the differencing stencil such that it accounts for the direction of wave prop- 
agation. Prior to the current investigation, there was one study into using this form 
of upwind differencing for the artificial compressibility method by Hartwich and Hsu 
[16,17]. However, they found it necessary to ensure that the scheme was total- variation 
diminishing (TVD) by reducing the order of the flux differencing and thereby adding 
more dissipation near sharp gradients. It was reasoned that since the incompressible 
Navier-Stokes equations will not produce solutions with strong discontinuities, such as 
shock waves, that flux- difference splitting could be applied without the use of TVD 
flux-limiters. 

In addition to the use of upwind differencing in the current formulation, an 
implicit solution algorithm was sought which would enhance the overall robustness 
of the code. The approximate factorization method of Beam and Warming [18] and 
Briley and MacDonald [19] has received much use since it’s introduction. However, 
the factorization error introduced can severely limit the time step. In addition, when 
applied to the method of artificial compressibility, it also limits the magnitude of 
the artificial compressibility constant [9,11]. In the current formulation the set of 
numerical equations are solved using a nonfactored line relaxation scheme similar to 
that employed by MacCormack [20]. This scheme can produce a solution close to 
the exact solution of a direct solver, without using much more time than a factored 
scheme. Its drawbacks are that it is not nearly as straight forward in its implementation 
as approximate factorization, and the line relaxation is inherently recursive, limiting 
the amount of optimization which can be done on vector processors. The current 
implementation trys to reduce these drawbacks by implementing the scheme in the 
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most efficient way possible. Since this research was being performed on a Cray 2 
computer, ample memory (256 million words of processor memory) was available. In 
all cases extra memory was used when it enabled a savings in computational time. 

In the following chapters, the details of the artificial compressibility scheme and 
its use in solving the incompressible Navier-Stokes equations for steady-state and time- 
dependent problems are given., The upwind-differencing scheme is detailed, and then 
the implicit solution procedure is described. Details of the boundary condition proce- 
dures are also included. Numerous computed results are presented with comparisons 
to analytical solutions and experimental results. The test problems include the flow 
in a driven cavity, the flow over a backward-facing step, the steady and unsteady flow 
over a circular cylinder, flow over an oscillating plate, flow through a one-dimensional 
inviscid channel with oscillating back pressure, the steady-state flow through a square 
duct with a 90° bend, and the flow through an artificial heart configuration with 
moving boundaries. 
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CHAPTER TWO 


GOVERNING EQUATIONS AND ARTIFICIAL COMPRESSIBILITY 


2.1 Introduction 


In this chapter the equations governing constant density, incompressible, viscous, 
Newtonian fluid flow are presented. The generalized coordinate transformation is 
presented, and then the artificial compressibility method is introduced. The details 
of its formulation for the solution of both steady-state and time-dependent problems 
are given. Finally, the numerical evaluations of the coordinate transformation metrics 
leading to a free-stream preserving system is detailed. The majority of the material 
presented here will be for the 3-D formulation only, as the 2-D system is an obvious 
subset of the 3-D system. The primary exceptions to this are the viscous fluxes and 
the eigensystem of the Jacobian of the convective flux vectors. The details of these 
two formulations for both 2-D and 3-D are given in Appendix A and Appendix B, 
respectively. 

2.2 Governing Equations 


The equations presented here have first been nondimensionalized using the fol- 
lowing quantities 


Ui 


Ui = 


— 


t = 


iu 


ref 


u ref 

P ~ Pref 

P U lef 


X r ef 


? T ij ~ 


P U le / 


X re f 
if — 


v 


— Re 


-1 


( 2 . 1 ) 


X re fU re f 

where Ui = u,v, or iu, the three Cartesian velocity components, for i = 1,2,3, respec- 
tively; where X{ = ®,y, or z, the three Cartesian spatial coordinates, for i = 1,2,3, 
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respectively; where p represents pressure, p represents density, t represents time, u is 
the kinematic viscosity, Re is the Reynolds number, Tij represents the viscous stress 
tensor, and the subscript ref denotes reference quantities. These reference quantities 
are generally chosen to be the freestream quantities for external flows, or the quantities 
at the the inflow boundary for an internal flow. The Navier-Stokes equations which 


govern incompressible, constant density flow are written in conservative form, with the 
tildes dropped for convenience, as 
du dv dw 

= 0 

( 2 . 2 ) 


where 
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For modeling turbulent flows, these equations represent Reynolds averaged quantities, 
and using the Boussinesq approximation for the Reynolds stress, the viscous stress 
tensor can be written in the form 


T «j 



(2.4) 


where u t is the turbulent eddy viscosity. As the primary thrust of this paper is con- 
cerned with algorithm development and testing, only laminar flow calculations will be 
considered and the determination of v t will not be discussed. However, the formulation 


9 



of the current algorithm will allow the use of a spatially and time varying viscosity. 
This will make it possible to implement an eddy viscosity turbulence model in the 
future. 

2.3 Coordinate Transformation 


The equations are transformed into generalized coordinates using 

t = t{x>V,z,t) 

r] = T](x,y,z,t) (2.5) 

C = C(®iSMiO 

in order to facilitate the use of time-varying, body-fitted coordinates for any specific 
type of geometry to be used in the numerical calculations. This results in the following 

d_ fU\ d_(V\ d_ (W\ =0 
m = ~ e ' ) ~ Tn (f ~ M ~ a< (s " = " r 

where r represents the right-hand side of the momentum equations, J is the Jacobian 


of the transformation and 
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The quantites U , V, and W are contravariant velocity components and the metrics of 
the coordinate transformation are represented here using 



The complexity of the differential form of the viscous fluxes varies considerably 
based on whether either of the following two assumptions can be made: that the flow 
computations are taking place on an orthogonal mesh; or that the viscosity is spatially 
constant (the flow is both laminar and assumed to be Newtonian). The viscous fluxes 
e vi fv> and g v based on any of the four combinations of the two above assumptions are 
presented in Appendix A. 

It is noted here that when using artificial compressibility to solve unsteady- 
flow problems on a time-varying mesh, it is important to perform this coordinate 
transformation before introducing the artificial compressibility relation, otherwise an 
incorrect right-hand side of the continuity equation will result in the divergence of 
velocity being nonzero. 

The following two sections deal with the formulation which makes it possible to 

t 

solve Eqs. (2.7) numerically using the method of artificial compressibility for either a 
steady-state or time-dependent problem. 


2.4 Steady-State Formulation 


The artificial compressibility relation is introduced by adding a time derivative 
of pressure to the continuity equation 


dp 

dr 


-/3V-« 


( 2 . 8 ) 
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In the steady-state formulation the equations are to be marched in a time-like fashion 
until the right-hand side f of Eq. (2.6) and the divergence of velocity approach zero. 
The time variable for this process no longer represents physical time and so in the 
momentum equations t is replaced with r, which can be thought of as a pseudo- time 
or iteration parameter. Combining Eq. (2.8) with the momentum equations gives the 
following system of equations 

f = - A(i -*)-£<*- *.) - £(* - 6.) = -k (2.9) 

where R is defined here as the residual vector of these equations and where 
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Since fast convergence is desired and accuracy in pseudo-time is not an issue, the 


pseudo- time derivative is replaced by an implicit Euler finite-difference formula, giving 



where superscript n denotes quantities at the nth pseudo-time iteration level. The 
right-hand side is linearized resulting in 

ts; /+ (!§) ](-p" + 1 -£") = -*" (2-n) 

where I is a 4 x 4 identity matrix and where D = JD. If this equation were solved 
exactly as it is, then for very large At this would become a Newton iteration for a 
steady-state solution. However, it is often not feasible to form the exact Jacobian of 
the residual vector R. Before these details are discussed, however, an equation similar 
to Eq. (2.11) for time-dependent problems will be developed. 

2.5 Time- Accurate Formulation 


In the time- accurate formulation the time derivatives in the momentum equations 
are differenced using a second-order, three-point, backward-difference formula 


3u n+1 _ 4 u n + u n - J 


_ _^ n +i 


where the superscript n denotes the quantities at time t = nAt and f is the right- 
hand side of the momentum equations given in Eq. (2.6). To solve Eq. (2.12) for a 
divergence free velocity at the n+1 time level, a pseudo-time level is introduced and is 
denoted by a superscript m. The equations are iteratively solved such that £ n+1,m+1 
approaches the new velocity u n+1 as the divergence of 'u n+1,m+1 approaches zero. To 
drive the divergence of this velocity to zero, the following artificial compressibility 
relation is introduced: 


At 


= -/3V .u n+1 ' m+1 


(2.13) 
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where r denotes pseudo-time and f3 is an artificial compressibility parameter. In this 
form it can be seen that the constants /? and At are not independent. However, they 
are kept separate here primarily because of the following reason. In the numerical 
equation which approximates this partial differential equation, (3 is moved inside the 
divergence operator, and the change in pressure becomes a nonlinear function of (3 
because of the use of upwind differencing. Therefore, in the numerical equation, f3 and 
At become independent. 

Combining Eq. (2.13) with the momentum equations gives 

It (D n+l 


= - J R n+1 - m+1 - — (1.5Z> n+1,m - 2 D n + 0.5-D” -1 ) 
A V 


(2.14) 


where D is the same vector defined in Eq. (2.10) and R is the same residual vector 
defined in Eq. (2.9). Also appearing in this equation is It r which is a diagonal matrix 
and J m which is a modified identity matrix given by 

I tT = diag 
J m = diag[0, 1, 1, 1] 

Finally, the residual term at the m+1 pseudo-time level is linearized giving the follow- 
ing equation in delta form 


1 1^> ^5 L5 

At’ At’ At’ At 


Ijr (dR\ 
6D J 


n-\-l ,m" 




(2.15) 


_ _^n+l,m _ 5 £)n+l,m _ + 0.5£> n_1 ) 

At 


As can be seen, Eq. (2.15) is very similar to the steady-state formulation given by 
Eq. (2.11). In a sense the time-accurate formulation requires the solution of a steady- 
state problem in order to advance one physical time step. Both systems of equations 
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will require the discretization of the same residual vector R. The derivatives of the 
viscous fluxes in this vector are approximated using second-order central differences. 
The formation of the convective fluxes is not such a simple matter and is the subject 
of the next chapter. 

2.6 Free-Stream Preservation 


The metrics of the coordinate transformation are evaluated from the following 
formulas (Ref. 21) 

* A 

Zx = y v z<; - y(Z v ( x = y^z v - y v Z£ 

?3/ = z v x < ~ z < x v Cy = z t x v ~ z v X( 

Cz = x vVC ~ x CVv Cz = x ( y v - x v y ( 

y* = y<H - y^z c it = ix x t + i y yt + izz t 

Vy = z < X ( -Z ( x c rjt = rj x x t + fj y y t + rj z z t 


Vz = x cV( - x ty< 
where i x = £ X /J, etc., and where 


A A A A 

Ct = (x x t + CyVt + CzZt 


dx 

*l=g~ v etc. 

Special care is taken in evaluating these metrics to ensure free-stream preserva- 
tion. When the grid is stationary this can be accomplished by evaluating the metrics 
using the following difference formula 

ix = ^{8r,y)nr,{S^z) - ih,(6(y)ft((S v z) 
where \i represents an average; i.e., 


— 0 ( a i+i,j,k + Oti-l,j,k) 
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and where 6 represents a central difference; i.e., 


= 2^ ai+l ^’ k ~ ai - 1 d.fc) 


Similar formulas are used for the other metric terms. 

When the grid is moving with time, an additional constraint must be placed on 
the evaluation of the metrics to ensure free-stream preservation. The geometric con- 
servation law of Thomas and Lombard [22] gives a formula for evaluating the Jacobian 
of the transformation at the n + 1 time level when using a first-order-implicit Euler 
scheme. A similar formula corresponding to the second-order, three-point, backward- 
difference scheme used in the current work is given by 


jn+l 


= 1.5 | y — 

\J n 


2 J 


n-1 


-At 


+ &v(v 0 + 


<<(<<)] } 


-i 


Using this formula to update the Jacobian will ensure free-stream preservation on a 
moving grid. This was tested on a number of grids, and it was found that it worked well. 
In particular, the moving grid for the artificial heart calculations presented in Chapter 
6 was found to produce residuals on the order of 10 -3 for freestream conditions when 
no correction for the grid velocities was made. However, when the above correction 
was made, the residuals were on the order of machine zero. 
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CHAPTER THREE 


UPWIND DIFFERENCING 
3.1 Introduction 

Upwind differencing is used in the present scheme as a means of following the 
propagation of the artificial waves introduced by the artificial compressibility. The 
upwind differencing thus provides a dissipative scheme which automatically suppresses 
any oscillations caused by the nonlinear convective fluxes. In addition, the upwind- 
differenced flux vector will contribute to terms on the main diagonal of the Jacobian 
of the residual, whereas a central-differenced flux vector would not. This will help 
make the implicit scheme nearly diagonally dominant and contribute greatly to the 
robustness of the code. Even though the upwind flux differences are more costly to 
form, the speed up in convergence does result in a significant savings in computing-time 
requirements. 

The upwind scheme is derived from one-dimensional (1-D) theory, and then is 
applied to each of the coordinate directions separately. Flux-difference splitting is used 
here to structure the differencing stencil based on the sign of the eigenvalues of the 
convective flux Jacobian. The scheme presented here was originally derived by Roe 
[12] as an approximate Riemann solver for the compressible gas dynamics equations. 

3.2 Flux-Difference Splitting 

The derivative of the convective flux in the £ direction is approximated by 

dE Ei+i/z — Ei-i/i 
dt * A £ 


(3.1) 
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where Ei+ 1/2 is a numerical flux and i is the discrete spatial index for the £ direction. 
The numerical flux is given by 


Ei+1/2 = g E{Di+\) + E(Di) - &+1/2] 


(3.2) 


where the <f>i +\/2 is a dissipation term. For 4>i+\/2 — 0 this represents a second-order 
central-difference scheme. A first-order upwind scheme is given by 


&+1/J = & E t+l/2 ~ AE i+l/2 


(3.3) 


where A E± is the flux difference across positive or negative traveling waves. The flux 
difference is computed as 


A Et +1/2 = A ± (D)AD i+1/2 


(3.4) 


where the A operator is given by 


AA+1/2 — Di+i — D{ 


The plus (minus) Jacobian matrix has only positive (negative) eigenvalues and is com- 
puted from 

A ± = XiAfX- 1 (3.5) 

where is composed of all the positive (negative) eigenvalues, where each individual 
positive (negative) eigenvalue is computed as proposed by Buning and Steger [23] 

Af = - ± ^A? + 

where e is a small parameter used to avoid a discontinuous transition in the delta 
fluxes as any of the eigenvalues changes sign. A value of e = 0.01 is used for all of the 
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calculations in this work. In Eq. (3.5) the subscript 1 denotes matrices corresponding 
to the ^-direction flux. The matrices X\ and X^ 1 are the right and left eigenvectors 
of the Jacobian matrix of the flux vector. These matrices for both the 2-D and 3- 
D systems of equations are presented in Appendix B. All matrices appearing in the 
upwind dissipation term must be evaluated at a half point (denoted by i+1/2). To 
do this a special averaging of the dependent variables at neighboring points must be 
performed. The Roe properties [12] which are necessary for a conservative scheme, are 
satisfied if the following averaging procedure is employed 

D = i(Z> i+ i + Di) (3.6) 

A scheme of arbitrary order may be derived using these flux differences. Im- 
plementation of higher-order-accurate approximations in an explicit scheme does not 
require significantly more computational time if the flux differences A are all com- 
puted at once for a single line. A third-order upwind flux difference can be obtained 
using 

&+1/2 = -g (& E t- 1/2 - ^ E t+ 1/2 + AE r+ 1/2 - AE r+ 3/2) ( 3 - 7 ) 

The primary problem with using schemes of accuracy greatei than third order 
occurs at the boundaries. Special treatment is needed, requiring a reduction of order 
or a much more complicated scheme. Therefore, when going to a higher-order-accurate 
scheme, compactness is desirable. Such a scheme was derived by Rai [24] using a fifth- 
order-accurate, upwind-biased stencil. A fifth-order, fully upwind difference would 
require 11 points, but this upwind-biased scheme requires only seven points. It is 
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given by 


*+ 1/2 = - 4 l_2AE * + -=/2 + 11 AB t ./2 - - 3 A *? + 3/2 

+2A£? {+5 ^ 2 — llAE i+3/2 + 6 A-E^ + j/ 2 + 

Next to the boundary, a scheme that is nearly second-order accurate can be 
maintaine 4 . by the third- and fifth-order schemes by using the following 

<t>i+l/2 = € (& E t+ 1/2 “ ( 3 ‘ 9 ) 

For € = 0, this flux leads to a second-order central difference. For e = 1, this is the 
same as the first-order dissipation term given by Eq. (3.3). By including a nonzero 
e, dissipation is added to the second-order, central-difference scheme to help suppress 
any oscillations. A value of e = 0.01 is used for all of the results presented in this 
paper. 

The right and left matrices given in Appendix B clearly show that the artifi- 
cial compressibility parameter /3 will affect not only the continuity equation, but the 
momentum equations as well. An analysis for the Cartesian coordinate case shows 
that the dissipation terms added to the momentum equations will grow as the square 
root of /?. This indicates that the value of (3 should be chosen with care when us- 
ing the upwind-differencing, as extremely large values could cause large errors in the 
differencing of the momentum equations. 

3.3 Central Differencing and Artificial Dissipation 

For some of the test cases to be presented in the Computed Results chapter a 
comparison between the present upwind-difference scheme and a central-differencing 
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plus artificial-dissipation formulation is made. This section presents the artificial terms 
used in these comparisons. 

As stated above, Eq. (3.2) will result in a central-difference scheme if = 0. 

If this formula is nondissipative, and if it is used without any additional numerical 
dissipation, an odd-even grid-point decoupling will develop, and the solution will show 
wiggles and oscillations from grid point to grid point. This is caused by the nonlinear 
convective fluxes introducing high frequencies which the natural viscous dissipation 
will not always damp out. A common method by which artificial dissipation can be 
added to the differencing is to use a difference approximation to a fourth derivative 
of the dependent variables on the right-hand side of the equations, and a difference 
approximation to a second derivative on the left-hand side. For more details, see 
Pulliam [25]. The right-hand side fourth-order dissipation can be given by 

0;+ 1/2 = — e e 1^4-1 |i+l/2(-^i-|-2 — 3-Dj-fi + 3 D{ — D{- 1) (3.10) 

where e e is the explicit artificial dissipation coefficient, and |A.i|; + i/ 2 is the spectral ra- 
dius of the Jacobian matrix of the convective flux vector. The second-order dissipation 
term used in forming the implicit side of the equations is given by 

4>i+ 1/2 = Ci|-Al|i + i/ 2 (-Dt+l — Di) (3-11) 

where e,- is the implicit artificial dissipation coefficient. 
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CHAPTER FOUR 


IMPLICIT SCHEME 

This chapter describes the way in which Eqs. (2.11) and (2.15) are numerically 
represented and solved for the 3-D equations. The 2-D system is solved in a similar 
fashion. The first consideration is the formation of the Jacobian matrix of the residual 
vector R required for the implicit side of the equation. Applying the upwind-difference 
formula given in Eq. (3.1) to the convective flux vectors, and applying a second- 
order, central- difference formula to the viscous terms, the residual at a discrete point 
(xijk,yijk,Zijk) is given by 


t3k A£ At ? 

~ 2 _ (-^'y)t+l/2,j,fc ~ (•^'v)i-l/3,j,fe 

A< A* 

(-^)»,j+l/2,fc ~ (-Fi>)i,j-l/2,fc _ (£t»)t,j,fc+l/2 ~ {Gv)i 1/2 

A77 A77 


(4.1) 


The generalized coordinates are chosen so that A£, A 77 , and A( are equal to one. To 
limit the band- width of the implicit system of equations the Jacobian of the residual 
vector will be formed by considering only first-order contributions to the upwind nu- 
merical fluxes as well as the second-order differencing of the viscous terms. Hence, the 
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cmly portion of the residual vector that is actually linearized is the following 


Rijk — 


2 (-^t+1 Ei -J,j,k + Fij +lt k — Fi t j-I t k + — Gi,j t k-1 

^ AF i+l/2,j,k + AF i-l/2,j,k ~ AF i-l/2,j,k 
~ AF ^/2,k + AF ^j+l/2,k + - 1/2>fc - ^-._ 1/2fJk 

A ^i,j,k+ 1/2 + A ^i,j,k+ 1/2 + A ^t,j,k- 1/2 ~ A ^i,j,k-l/2) 


~(Fv)i+l/2,j,k + (#i>)t -1/2, j , k ~ ( F v )i,j+l/2,k 

+(F V )iJ-l/2,k ~ {G v )i,j,k+l/2 + {G v )i,j,k- 1/2 


(4.2) 


The exact Jacobian of this residual vector will result in a banded matrix of the form: 


dR 

8D 




dRijk 0 0 - dRijk 0 0 
, U, U, n ^ , u, u, 


[dD iijtk -i ’ ’ 'dDij-u 

dRjjk dRjjk dRjjk 
dDi-i,j,k &Dij } k dDi+ .'ij'k 
0 0 _9Rjjk _ _ Q . Q dR ijk 

U ? —> U > OT-k 5 U > — J 






(4.3) 


where each entry of the banded matrix represents a vector of 4 x 4 blocks in 3-D (3x3 
blocks in 2-D) which is aligned along a diagonal of the matrix. These exact Jacobians, 
however, can be very costly to form. Therefore, approximate Jacobians of the flux 
differences as derived, and analyzed by both Barth [26] and Yee[27] are used. These 
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are given as follows 


dDi j k-l ^t,j , fc — 1/2 + ^i,j, fc-1/2) (T3 )i,7,* — 1/2 

dD~ 1)fc ~2 ~ ^.1-1/2,* + B iJ~l/2,k) ~ Mi,j-l/2,k 

dD~^~k ~ ^t-i/2 ,j,k + A-i/2,j,fc) " (71)1-1/2,7,* 

4 + 4- 4+ _ ,4 - _ 4 — 

dDij'k 2^ *+1/2,7,*^ n i-\/2,j,k A t+l/2,7,fc A i-i/2j,k 
+ -®ij+l/2,fc + B t,j-l/2,k ~ A, 7+1/2,* - A, 7-1/2,* 

+ *+l/2 + fc-1/2 ~~ A,j,fc+ 1/2 ~ ^t.j, fc-l/2) 

+ (“Tl )i+l/2,j,fc + {j2)i,j+l/2,k + (73 )i,j,fc+l /2 
+ (71)1-1/2,7,* + (72)i,j-l/2,* + (73)i, 7, *-1/2 
'^(-A+l,/,* — A+l/2,7,fc A+l/2,7,fc) — (7l)i+l/2,7,* 


dR 


ijk 


dDi 


*+iii>fc 


dRijk 1 . - + _ 

~o ‘^‘’• 7+1,fc ~ A, 7+1/2,* + A, 7 + 1/2,*) _ ( 72 )i, 7+1/2,* 


(4.4) 


0A.7+1,* 2 

0A,J,*+1 ~ 2^ i,i,fc+1 ~ ^£.*+1/2 + ^iTy.fc+i/z) - (73)1,7’, *+1/2 
where ^4 = Aii-B = A2, and C = A3 as given in Appendix B, and only the orthogonal 
mesh terms are retained for the implicit viscous terms giving 

7i 


1 (£ + ^+£)/mjj 


ReJ 


72 


1 (vl+vl+vl)!™-^ 


ReJ 


73 ReJ + A + A )■ Im 

(Tl )i+l/2,j,A: ^ [(Tl )i+l "I" (Tl )i — 3 


(4.5) 


(72)i, 7+1/2,* - 2 [(72)i,7 + l,* + (72)1,7-1,*] 


(73)t,7,fe+l/2 - 2 [(73 )»,7,*+l + (73)i,7,*-l] 
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The term Im is a modified identity matrix with a zero in the first diagonal entry. 

The matrix equation is solved using a Gauss-Seidel line- relaxation method similar 
to that used by MacCormack [20] and by Chakravarthy[28]. To implement this scheme, 
the entire numerical matrix equation is first formed from values at the nth time level. 
This takes a significant amount of memory to store all of the left-hand side matrices, 
but the savings in computing time justifies this. At this point the numerical equation 
is stored as the following banded matrix 


B[T, 0, ..., 0, U, 0, ..., 0, X , F, Z, 0, ..., 0, V, 0, ..., 0, W]AD = R 


(4.6) 


where AD ~ D nJrl — D n and T, U, V 9 W, X, Y, and Z are vectors of 4 x 4 blocks 
which lie on the diagonals of the banded matrix, with the Y vector on the main 
diagonal. This matrix equation is approximately solved using an iterative approach. 
One family of lines is used as the sweep direction. Using, for example, the £ family, a 
tridiagonal matrix is formed by multiplying the elements outside the tridiagonal band 
by the current AD and shifting them over to the right-hand side. This system can be 
represented by the following for a forward sweep 

BIX,Y,Z)&D‘» = R~TAD‘"_ 1 - UAD'+'_ hk 

~ VAD‘ J+lti - \VAD\jw 

where l is an iteration index. The number of iterations is generally a small number like 
two or three. This system can be solved most efficiently by first performing and storing 
the lower-upper (LU) decomposition of the tridiagonal matrix before the iteration is 
begun. Then for each iteration, the right-hand side is formed using the latest known 
A D (which is set to zero for the first iteration), and the entire system is backsolved. 
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The LU decomposition can be entirely vectorized, but the back solution is inherently 
recursive and cannot be vectorized. 
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CHAPTER FIVE 


BOUNDARY CONDITIONS 

5.1 Method of Characteristics 


Implicit boundary conditions arc used at all of the boundaries, this helps make 
possible the use of large time steps. At a viscous no-slip surface, the velocity is 
specified to be zero, and the pressure at the boundary is obtained by specifying that 
the pressure gradient normal to the wall be zero. The boundary conditions used for 
inflow and outflow regions are based on the method of characteristics. The formulation 
of these boundary conditions is similar to that given by Merkle and Tsai [ 29 ], but the 
implementation is slightly different. The scheme is derived here for a constant £ 
boundary, with similar results for a constant 77 or a constant £ boundary. The finite- 
speed waves which arise with the use of artificial compressibility are governed by the 
following 


dD _ dE _ dEdD _ -dD J0D 

dr ~ at ~ dD dt ~~ A dz ~ ~ XAX ajf 

Multiplying by A -1 gives 


x ~'t = - AX ”f 


( 5 . 1 ) 


If one were to move the X -1 matrix inside the spatial and time derivatives, the re- 
sult would be a system of independent scalar equations, each having the form of a 
wave equation. These equations are termed the characteristic equations. The sign 
of the eigenvalues in the A matrix determines the direction of travel of each of the 
characteristic waves. For a positive (negative) eigenvalue, there corresponds a wave 
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propagating information in the positive (negative) £ direction. The number of positive 
or negative eigenvalues determines the number of characteristic equations propagating 
information from the interior of the computational domain to the boundary. Thus, 
at the boundary, we will use these characteristics which bring information from the 
interior as part of our boundary conditions. The rest of the information should come 
from outside the computational domain, and so one is free to specify some variables. 

There will always be at least one characteristic wave traveling toward the bound- 
ary from the interior because there is always at least one positive eigenvalue and one 
negative eigenvalue in the artificial compressibility equations. To select the proper 
characteristic waves, Eq. (5.1) is multiplied by a diagonal selection matrix L which 
has an entry of one in the position of the eigenvalue we wish to select, and zeros 
elsewhere. Thus 


LX ~'f = - £Ajr_, f 


(5.2) 


Replacing the time derivative with an implicit Euler time step gives 


( W + £AX "‘ I) {D " + ' ~ Dn) = ~ LAX ~' d ~W (5 - 3) 

The following discussion relates to the number of positive and negative eigen- 
values at each specific type of boundary for the 3-D system of equations. Similar 
results apply to the 2-D system, the only difference being that where there are three 
eigenvalues of one sign in the 3-D system, there are only two eigenvalues in the 2-D 
system. 

Equation (5.3) gives either one or three relations, depending on the number of 
nonzero elements in L. To complete the set of equations, some other information must 
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be supplied. Here is defined a vector Cl that contains a boundary condition for each 
zero element on the diagonal of the matrix L. The other elements of Cl are set to zero. 
Since these boundary conditions are held constant in time, the following holds true 


dCl n 
dr ~ 


dCl dD 
8D &t 


= 0 


dCl 

dD 


= 0 


(5.4) 


Combining Eqs. (5.3) and (5.4) gives 

LX - 1 


+ LAX - 1 4 + ^1 (D n+1 - D n ) = -LAX- 1 ^ 


(5.5) 


JAr 1 d£ dD J 

Equation (5.5) can be used to update the variables implicitly at any of the inflow or 
outflow boundaries with the proper choice of L and f). 


5.2 Inflow Boundary 


At the inflow, there will be one characteristic wave traveling out of the compu- 
tational domain and three traveling in since fluid is traveling into the domain. If the 
incoming fluid is traveling in the positive £ direction, then the signs of the eigenvalues 
are given by 

A l5 A 2 >0 

A3 > 0 

A 4 < 0 

This fourth eigenvalue corresponds to the only wave carrying information out of the 
computational domain. So L will have a one for the fourth diagonal entry. If the 
incoming fluid is traveling in the negative £ direction, then 

Ai, A 2 <0 

A3 > 0 
A 4 < 0 
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For the same reasoning a one in the third diagonal entry of L is required for this 
boundary. 

Two different sets of specified variables have been used successfully for inflow 
boundaries. One set consists of the total pressure and the cross-flow velocity. This set 
is useful for problems in which the inflow velocity profile is not known. For this set 
and for the first boundary mentioned previously, the ft vector is 

1 U V w~ 

0 0 10 
0 0 0 1 
0 0 0 0 . 

The second possible set of specified variables consists of the velocity components. This 
is useful for problems in which a specific velocity profile is desired at the inflow. The 
ft vector for this is 

oioo- 
0 0 10 
0 0 0 1 
0 0 0 0 . 

5.3 Outflow Boundary 

At the outflow boundary there are three characteristic waves traveling out of the 
computational domain and one traveling in since fluid is also leaving the domain. If 
the fluid is traveling along the positive £ direction then 

Aj , A 2 >0 
A3 > 0 
A4 < 0 

and ones are required in the first three diagonal entries of the L matrix. If the fluid is 


■ u 


V 

^ft 

w 

5 dD ~ 

.0. 



"p + +V + ™ ) 


V 

5 ft 

w 

; dD ~ 

0 
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traveling in the negative £ direction then 


Ai , A 2 <0 

A3 > 0 

A 4 < 0 

and ones are required in the first, second, and fourth diagonal entries of the L matrix. 

For all of the test problems presented in this paper static pressure was specified 
at the outflow boundary. For the first outflow boundary mentioned previously this 
results in 

-0000- 
0 0 0 0 
0 0 0 0 
.1 0 0 0 . 
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CHAPTER SIX 


COMPUTED RESULTS 
6.1 Introduction 

The 2-D and 3-D codes have been run for a large number of laminar flow test 
cases. Eight of these cases are presented here as examples of the codes capability to 
efficiently and accurately predict flow physics. These cases are as follows: a driven 
square cavity flow, which has become a standard validation test case due to its simple 
boundary conditions; the flow .over a backward facing step, which has also become 
a popular validation case as it is an example of an internal flow with a recirculating 
region; the steady flow over a circular cylinder, which has become an extensively 
validated external flow problem; the flow on top of an oscillating plate, known as 
Stoke’s second problem; the flow through an inviscid one dimensional (1-D) channel 
with an oscillating back pressure; vortex shedding behind a 2-D circular cylinder; the 
flow through a 3-D square duct with a 90° bend; and the flow through an artificial 
heart. 

The computing times reported here are the CPU seconds used on a Cray 2. For 
comparison, these times are nearly the same as obtained running on a Cray XMP-48. 
For steady-state problems, the computations are run until the maximum residual has 
converged between 4 to 6 orders of magnitude, the maximum divergence of velocity over 
all the points is less than 10 -4 , and the flow quantities being measured have approached 
a steady-state value in at least 4 significant digits. For unsteady-flow problems the sub- 
iterations were continued until the maximum residual and the maximum divergence of 
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velocity was less than 10 -4 , except for the artificial heart case, in which this number 
was increased to 10 -3 due to computing time limitations. This convergence is usually 
obtained in 100 to 200 iterations for steady-state cases, and in 10 to 20 subiterations 
for unsteady problems. 


For each of the test cases presented, the larger the time step At, the better 
the convergence was. In all of the upwind-difference cases presented here the solution 
remained stable no matter how large a time step was used, so the time step was set to 
10 12 which effectively reduced the 1/Ar term to zero. The choice of f3 for each case was 
determined through numerical convergence tests. It was found that the convergence 
was quite sensitive to the value of /?, and in some cases, the choice of /3 could cause 
the scheme to become unstable. For most cases, however, it was easy to find a range 
of /5 for which the code would converge very quickly. This process usually required 
running the code for a small number of iterations using values of (3 equal to 0.1, 1, 
10, 100, etc., and isolating the range of (3 which gave the best convergence rates. The 
convergence of the current formulation is degraded by the errors introduced by the 
approximate Jacobians on the left-hand side of the equations and by the fact that the 
whole system of equations is not exactly solved by the line-relaxation process. If it 
were possible to use the exact Jacobians and solve the system exactly, then this would 
be a Newton iteration, in which case one would expect to have quadratic convergence 
when using a very large time step for any value of /?. In all cases, the time quantities 
given in the following flow problems refer to dimensionless time. 
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6.2 Driven Cavity Flow 


The two-dimensional flow in a driven square cavity whose top wall moves with a 
uniform velocity has been used rather extensively as a validation test case by several 
authors in the recent past. It provides a good test case in that there is no primary 
flow direction and the boundary conditions are very simple to employ. Ghia et al. 

[30] presented extensive numerical data obtained from their multigrid vorticity-stream 
function formulation using very fine grids. Although there was no experimental results 
for this problem, they reported results which agreed quite well with other computa- 
tional efforts, particularly at the lower range of Reynolds numbers. Other recent 
computational work involving this particular geometry include Schreiber and Keller 

[31] who use a vorticity-stream function formulation; Kim and Moin [6] who use a 
fractional-step method in primitive variables in conjunction with approximate factor- 
ization; Vanka [32] who uses a multigrid technique in primitive variables; and Benjamin 
and Denny [33] who use a centrally-differenced vorticity-stream function formulation 
in conjunction with an ADI scheme. 

The current calculations attempt to maintain the accuracy of these authors while 
using fewer grid points. First, the results of a grid resolution study are reported here. 
This study was done for a Reynolds number of 10,000, which was the highest Reynolds 
number tested and had by far the most complex flow pattern. The grids for this study 
were obtained by first computing a grid with 161 by 161 points, which were clustered 
near the walls. Subsequently coarser grids of dimensions 81 by 81, 41 by 41, and 21 
by 21 were obtained by removing every other grid point in both directions from the 
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previous grid. An example of one of these grids, the 81 by 81 grid, is shown in Fig. 1. 
These calculations used 11 implicit line relaxation sweeps in the ^-direction for each 
iteration. The parameters used in these calculations, as well as several flow quantities 
are shown in Table 1. The value of /?, the number of iterations required for 6 orders 
of magnitude convergence ( nt c ), and the computing time, comprise the first entries in 
this table. This shows how the value of the near optimum /? for this problem increases 
with the grid-point density. The convergence does slow somewhat with increasing 
grid-point density, but it is only by a small amount for the denser grids. The other 
entries in Table 1 are several of the flow quantities which were computed, namely the 
maximum and minimum vertical velocity component along the centerline of the cavity, 
and the streamfunction, V’mtnj and vorticity values, u> v . c ., at the center of the primary 
vortex. 


Table 1 Numerical Parameters and Flow Quantities for 
Grid- Resolution Study of Driven-Cavity Flow 


Grid 

/? 

nt c 

Time 

%fll 


i’min 

W».c. 

161 

1 

440 

680 

0.458 

-0.573 

-0.1204 

-1.880 

81 

1 

390 

154 

0.446 

-0.562 

-0.1174 

-1.817 

41 

0.1 

380 

42 

0.370. 

-0.473 

-0.1015 

-1.621 

21 

0.01 

170 

6 

0.059 

-0.070 

-”.0222 

-0.790 


Table 2 shows the relative error in each of these flow quantities for the coarser 
grids as compared to the 161 by 161 grid results. Also included in this table are 
the ratios of the error from one grid to the next. These values show the order of the 
convergence toward a grid-independent solution. Since these values are all greater than 
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Fig. 1. Grid with 81x81 points used for computing the driven cavity flow. 

4, this indicates that accuracy of the overall spatial differencing is between second and 
third order. Thus the effect of the fifth-order upwind differencing is to help increase 
the overall order above the second-order accuracy of the viscous terms and the near- 
boundary differencing. These numbers indicate that the difference between the 161 by 
161 grid and a finer grid would be very slight, and that the solution given by this grid 
is essentially grid independent. This also shows that the results from the 81 by 81 grid 
are only in error by 2 or 3 percent, and that the results from this grid are close enough 
to a grid-independent solution that this grid could be used for the rest of this study. 

The remaining flow calculations for this problem were computed using the 81 by 
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Table 2 Relative Errors in Flow Quantities for 
Grid-Resolution Study of Driven-Cavity Flow 


Grid 

Vina x 

v min 

^mtn 

Wv.c, 

81 

0.024 

0.019 

0.025 

0.034 

41 

0.191 

0.173 

0.157 

0.138 

21 

0.872 

0.878 

0.816 

0.580 

er4i/er 8 i 

7.96 

9.10 

6.30 

4.06 

er 2 i/er 41 

4.57 

5.08 

5.20 

4.20 


81 grid for Reynolds numbers of 100, 400, 1000, 3200, 5000, 7500, and 10,000. The 
value of the artificial compressibility (3 was set to 20 for the Reynolds number of 100, 
to 10 for the 400 Reynolds number case, to 2 for a Reynolds number of 1000, and 
was set to 1 for the higher Reynolds numbers. This problem indicates that there is an 
inverse relationship between the Reynolds number and the optimum (3. 

The velocity components on the fines passing through the geometric center of 
the cavity are compared to the results of Ghia et al. [30] in Fig. 2. The u-velocity 
component is plotted along the y-axis for the different Reynolds numbers in Fig. 2a. 
The origins of the plots has been shifted to the left for each successive Reynolds number 
case. The data of Ghia was obtained from a uniform grid of 129 by 129 points for 
Reynolds numbers up to 3200, and a uniform grid of 257 by 257 points for Reynolds 
numbers 5000 and greater. It is seen that these two computed results agree very well 
with each other. In Fig. 2b, the v-velocity component is plotted along the x-axis 
passing through the geometric center for the different Reynolds numbers. The origins 
of these plots are shifted up for each successive Reynolds number case. Again, good 
agreement is seen between the two computed results. 
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2a. U- velocity component versus y on the vertical centerline. 



X 


2b. V-velocity component versus x on the horizontal centerline. 

Fig. 2. Comparison between present results (solid line) and computations of Ghia et 
aZ. [30] (symbols). □ : Re=100, O* Re=400, A: Re=1000, +: Re=3200, X: Re=5000, 
<C>: Re=7500, and y: Re=10,000. 


In Table 3, the stream function and vorticity quantities are given for the core 
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Table 3 Stream-Function and Vorticity at the Center of the 
Primary Vortex for Different Reynolds Numbers 


Re 

Present 
V’Tnin (^v.c.) 

Ghia et al [24] 

^mtn K.C.) 

Schreiber et al. [25] 

V’mtn ( w t>.c.) 

Kim et al [6] 

rffnin ( w #.c.) 

100 

-0.1030(-3.104) 

81x81 

-0.1034(-3.166) 

129x129 

-0.1033(-3.182) 

121x121 

-0.1030(-3.177) 

65x65 

400 

-0.1131(-2.296) 

81x81 

-0.1139(-2.294) 

129x129 

-0.1130(-2.281) 

141x141 

-0.1120(-2.260) 

65x65 

1000 

-0.1171(-2.044) 

81x81 

-0.1179(-2.050) 

129x129 

-0.1160(-2.026) 

141x141 

-0.1160(-2.026) 

97x97 

3200 

-0.1195(-1.904) 

81x81 

-0.1204(-1.989) 

129x129 

- 

-0.1150(-1.901) 

97x97 

5000 

-0.1192(-1.846) 

81x81 

-0.1190(-1.860) 

257x257 

- 

-0.1120(-1.812) 

97x97 

7500 

-0.1186(-1.846) 

81x81 

-0.1200(-1.880) 

257x257 

- 

- 

10000 

-0.1177(-1.826) 

81x81 

-0.1197(-1.881) 

257x257 

-0.1028(-1.622) 

180x180 

— 


of the primary vortex for all Reynolds numbers tested. Included with the present 


results are the results of Ghia et al [30], Schreiber and Keller [31], and Kim and 
Moin [6]. Listed below the flow quantities is the grid size used for the calculation. 
Good agreement among all calculations is seen in the lower Reynolds number cases. 
The discrepancies between different solutions increase at the higher Reynolds numbers, 
although the same general trend of a leveling off and then a slight decrease in the value 
of the stream-function is seen. 

To study the 10,000 Reynolds number case in more detail, the streamlines are 
plotted in Fig. 3. The values of the stream-function contours for this plot are given 
in Table 4. The contour levels plotted correspond with the values plotted by Ghia 
et al. [30] for this case. Qualitatively, the plots appear to be identical. They each 
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Contour 

number 

Value of ip 

Contour 

letter 

Value of ip 

3 

1.0 x 10“ 4 

A 

-1.0 x 10- 5 

4 

2.5 x 10~ 4 

B 

-1.0 x 10“ 4 

5 

5.0 x 10 -4 

C 

-0.01 

6 

1.0 x 10“ 3 

D 

-0.03 

7 

1.5 x 10~ 3 

E 

-0.05 

8 

3.0 x 10“ 3 

F 

-0.07 



G 

-0.09 



H 

-0.1 



I 

-0.11 



J 

-0.115 


show secondary vortices of the same size and shape in the lower corners and the upper 
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left corner. In Table 5, the stream-function V’v.c.j vorticity w». c> , and location of the 
vortex core (x v , c .,y v .e.) for all the secondary vortices for this 10,000 Reynolds number 
case are given for both the present results and the results of Ghia et al. [30]. In this 
table, the initial T stands for top, B stands for bottom, R stands for right, L stands 
for left, and the superscript number corresponds to the level of the secondary vortex. 
Thus BR 3 refers to the third and smallest secondary vortex found in the bottom right 
corner of the cavity, which is not seen in Pig. 3 because it is too small. Quite good 
agreement between the two computations is seen for this case, especially considering 
that the results of Ghia et al. [30] use over 10 times as many grid points (66,049 versus 
6561). 


Table 5 Properties of the Secondary Vortices for the 
Driven Cavity at Re = 10,000 


Vortex 

Results 

ifiv.c. 

^V.C. 


Yv.c. 

TL 

Present 

2.418xl0 -3 

2.191 

0.0723 

0.9117 


Ghia et al. [30] 

2.420x1 0~ 3 

2.183 

0.0703 

0.9141 

BL 


1.434xl0~ 3 

2.097 

0.0585 

0.1686 



1.518xl0 -3 

2.086 

0.0586 

0.1641 

BR 


3.227x1 0~ 3 

4.163 

0.7619 

0.0585 



3.418xl0~ 3 

4.053 

0.7656 

0.0586 

BL 2 


-5.120xl0“ 7 

-0.02207 

0.1416 

0.0172 



-7.757xl0" 7 

-0.02754 

01560 

0.0195 

BR 2 


-2.103xl0 -4 

0.3726 

0.9277 

0.0729 



-1.313xl0 -4 

0.3126 

0.9336 

0.0625 

BR 3 


-4.267x10 “ 7 

-2.956xl0 -3 

0.9981 

0.0087 



-5.668xl0~ 9 

- 

0.9961 

0.0039 


The convergence toward a steady-state for this problem was very good for the 
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3 lowest Reynolds number cases, which required less than 150 iterations and only 
50 seconds of computing time. The higher Reynolds number cases were slower to 
converge, the 10,000 case took 390 iterations and 150 seconds of computing time. The 
average of the computing requirements for all seven cases came out to 200 iterations 
and 75 seconds of computing time. 

6.3 Plow Over a Backward Facing Step 

A second two-dimensional problem which has been used as a validation case is 
the flow over a backward-facing step. This case was run using both upwind differenc- 
ing and central differencing. The challenge in modeling this problem comes from the 
fact that the size of the separation bubbles downstream of the step are very sensitive 
to the pressure gradient in the flow. The geometry used in these calculations is shown 
in Fig. 4. At the inflow boundary, a parabolic profile is prescribed throughout the 
calculation, and the static pressure is allowed to change. Two step-heights downstream 
from the inflow a one-to-two expansion is encountered. The outflow boundary extends 
to 30 step heights downstream of the step, where characteristic boundary conditions 
are used along with the specification that the static pressure remain constant. The 
ability of the flow code to predict the reattachment length, xl, of the primary separa- 
tion bubble behind the step as well as the separation and reattachment locations, x2 
and x3, of the secondary separation bubble on the opposite wall was tested by com- 
paring the computed results from both the upwind- and central-differencing schemes 
to experimental values given by Armaly et at. [34]. 
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6.3.1 Upwind-Differencing Calculations 

The upwind difference scheme was used to compute these quantities for the 
laminar range of Reynolds numbers from 100 to 800, which are based on the average 
inflow velocity and twice the step height. The flow was calculated using a grid of 100 
points in the streamwise direction and 53 points in the cross-flow direction downstream 
of the step. The streamwise points were clustered toward the vertical-step face, and 
the cross-flow points were clustered at the walls. The initial conditions were specified 
to be freestream velocity at the interior points with uniform pressure everywhere. For 
the Reynolds numbers of 100 and 200, (3 was set to 1, for the Reynolds number of 300 
case, 0.5 was used, and for the Reynolds numbers of 400 through 800, / 3 was set to 0.1. 
The implicit line-relaxation process used 2 sweeps along the primary flow direction. 

In Fig. 5, the quantities xl, x2, and x3 are plotted versus Reynolds number for 
both the present computed results and the experimental results of Armaly et al. [34]. 
Good agreement is seen between the two for the value of xl at the lower Reynolds 
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01 




200 400 600 800 

Reynolds Number 


100 


Fig. 5. Separation length versus Reynolds number for the flow over a backward facing 
step. Solid line: computed xl, dashed line: computed x2, dotted line: computed x3, 
A: experimental xl, o : experimental x2, and o : experimental x3. 


numbers before the secondary separation appears. At a Reynolds number of 400, 
the secondary separation bubble is present, and the computed primary reattachment 
length begins to fall off of the experimental curve. Similarly, the computed secondary 
separation distance is shorter than the experimental values, although they do agree in 
that the secondary separation point is upstream of the primary reattachment point. 
The computed secondary reattachment point is seen to be close the experimental 
values. In their experiment, Armaly et al. [34] reported that the flow was found 
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to be three-dimensional near the step for Reynolds numbers greater than 400, and 
that the three-dimensional effects were negligible for lower Reynolds numbers. These 
three-dimensional effects could explain the discrepancies between the calculations and 
experiment. 

Results similar to the present results were reported by Kim and Moin [6]. They 
reported a primary reattachment length of just under 12 step-heights for a Reynolds 
number of 800, and the present result for this Reynolds number is 11.48. They reported 
a secondary separation bubble size of 7.8 and 11.5 step-heights for Reynolds numbers 
of 600 and 800, respectively. The present results show secondary separation bubble 
sizes of 7.34 and 11.07 step-heights for these two Reynolds numbers. The similarities 
between these computational results give additional credence to the idea that the 
three-dimensionality of the experiment of Armaly et al. [34] was the reason for most 
of the discrepancies exhibited in Fig. 5. 

The convergence characteristics of the upwind code for this problem are very 
good. In Fig. 6 the convergence histories of the Re = 100 and 800 cases are plotted. 
Figure 6a plots the log of the maximum residual normalized by the maximum residual 
at the first time step versus iteration number for the Re = 100 and 800 cases. Figure 6b 
plots the primary reattachment length xl versus iteration number. The Re=100 case 
converges within 55 iterations and the Re=800 case converges within 165 iterations. 
The average number of iterations required for the 8 different Reynolds number cases 
is 104 and the average required computing time is just under 11.5 seconds. 
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6.3.2 Central- Difference Calculations 


The central- difference scheme was used to compute the flow over a backward- 
facing step for a Reynolds number of 100, using the same grid, and the same initial 
and boundary conditions as the upwind-difference calculations. The quantity /3 was 
set to 2, and the pseudo-time step At was set to 0.11. This was the largest value 
of At for which the code would remain stable. The artificial dissipation coefficients 
e e and e; were set to 0.04 and 0.12, respectively. The separation bubble length was 
computed to be 2.84, which is less than 2% off the upwind- calculation value of 2.89. 
In general, the computational results from both methods agreed very well. 

Although the final solution was nearly the same as the upwind case, it took 1675 
iterations to converge, which took 164 seconds of computing time. So although the 
central difference flux vectors save about 15% of the total computing time cost per 
iteration, the overall computing costs are much greater. In addition, running this case 
took significantly more time by the investigator because of the need to specify the 
smoothing coefficients, as well as the maximum allowable time step, through trial and 
error testing. Since central differencing leads to equations which are far less diagonally 
dominant, the system is far less robust.' It also means that a line-relaxation scheme 
may not be the best way to solve the system of equations. 

Previous work by the author [11] in computing flow through this same geometry 
with the INS3D code gave similar convergence results as those seen here, and slightly 
better computing time requirements. The INS3D code uses central differences and 
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approximate factorization, which indicates that approximate factorization may be a 
more appropriate solution scheme in conjunction with central differencing. ' 

6.4 Steady Flow Over a Circular Cylinder 

As an example of an external flow problem, the flow over a two-dimensional 
circular cylinder was calculated using both upwind and central differencing. The grid 
was an algebraically generated O-grid with 100 points in the circumferential direction 
and 60 points radially. The grid points were clustered in the radial direction toward 
the body, and the outer boundary was placed 20 diameters from the cylinder. The 
code was run and steady-state solutions were obtained for Reynolds numbers of 5, 10, 
20, and 40, based on the freest ream velocity and the cylinder diameter. At the outer 
boundary method of characteristic boundary conditions were employed. In addition, 
where the fluid was entering the domain, the velocity was held constant, and where 
fluid was leaving the domain, the static pressure was held constant. 

For each case, several flow quantities were computed. Fig. 7 shows a schematic 
diagram of the geometry for this flow problem along with several of these flow quan- 
tities. These quantities are the flow separation length measured from the rear of the 
cylinder in cylinder diameters (i se p)> the angle which defines the point of separation 
from the body ( 6 3ep ), the coefficient of drag (Cp), the coefficient of pressure drag 
(Cd p ), and the coefficient of pressure at the front ( C p / ) and rear (C pr ) stagnation 
points. 

6.4.1 Upwind-Differencing Calculations 

In using the upwind-differencing scheme, the value of the artificial compressibil- 


48 



Fig. 7. Schematic diagram showing flow quantities for the circular cylinder flow 
computations. 


ity constant /? was set to 50 for all the cases. The line relaxation scheme used 4 sweeps 
in both coordinate directions, which seemed to be the best trade-off of convergence 
versus computing time. In Table 6, the flow quantities are presented for the present 
calculations as well as the numerical results of Takami and Keller [35], Dennis and 
Chang [36], Tuann and Olson [37], Braza et al. [38], and the experimental results of 
Coutanceau and Bouard [39] and of Tritton [40]. The comparisons show that there 
is very good agreement among nearly all the results. The results of the experiment, 
however, do not agree well with the computations. This could be because it is difficult 
to obtain accurate experimental results at such low Reynolds numbers. The conver- 
gence of the code toward a steady-state solution for this case was found to be quite 
good. All four Reynolds number cases converged in less than 70 iterations, requiring 
an average of 21 seconds of computing time. 

6.4.2 Central- Differencing Calculations 

The Reynolds number 5 case was run using the central-differencing scheme. The 
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Table 6 Flow Quantities for a Circular Cylinder 

Reynolds Number 


Source 

5 10 20 40 


L»ep 


Present 

0 

0.254 

0.932 

2.29 

Ref. 35 

0 

0.249 

0.935 

2.32 

Ref. 36 

0 

0.252 

0.94 

2.35 

Ref. 37 

0 

0.25 

0.9 

2.1 

Ref. 39 (exp) 

- 

0.34 

0.93 

2.13 



0,e P (degrees) 


Present 

0 

28.8 

43.1 

53.0 

Ref. 35 

— 

29.3 

43.7 

53.6 

Ref. 36 

0 

29.6 

43.7 

53.8 

Ref. 37 

>6 

29.7 

44.1 

54.8 

Ref. 38 

— 

- 

43.6 

54.5 

Ref. 39 (exp) 

- 

32.5 

44.8 

53.5 



Cd (Cd p ) 


Present 

4.18 (2.19) 

2.89 (1.602) 

2.08 (1.242) 

1.549 (1.011) 

Ref. 35 

— 

2.80 

2.01 

1.536 

Ref. 36 

4.12 (2.20) 

2.85 (1.600) 

2.05 (1.233) 

1.522 (0.998) 

Ref. 37 

4.66 (2.48) 

3.18 (1.775) 

2.25 (1.35 ) 

1.675 (1.095) 

Ref. 38 

— 

- 

2.18 

1.60 

Ref. 40 (exp) 

4.16 

3.06 

2.02 

1.65 



Gpf (- 

~Cpr) 


Present 

1.847 (1.067) 

1.476 (0.755) 

1.265 (0.615) 

1.147 (0.536) 

Ref. 35 

— 

1.474 (0.670) 

1.261 (0.537) 

1.141 (0.512) 

Ref. 36 

1.872 (1.044) 

1.489 (0.742) 

1.269 (0.589) 

1.144 (0.509) 

Ref. 37 

2.23 (1.081) 

1.744 (0.773) 

1.457 (0.614) 

1.312 (0.543) 


artificial dissipation coefficients c e and £; were set to 0.01 and 0.03, respectively. The 
line-relaxation scheme used 2 sweeps in the circumferential direction, and the artificial 
compressibility constant /3 was set to 50. The largest pseudo-time step At for which 
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the code remained stable was 0.01. This lead to very slow convergence. After 2000 
iterations the solution was very nearly converged. Again, the solution was very close 
to that of the upwind- differencing calculations. The computing time was about 260 
seconds, close to 10 times that required by the upwind scheme for this case. 

6.6 Oscillating Plate 

As an initial unsteady test case, the motion of fluid on top of an infinite plate 
which oscillates back and forth its own plane was calculated. This problem, also known 
as Stoke’s second problem (see for example White [41]), was set up with the x-axis 
along the plate, and the y-axis normal to it. The velocity of the plate is set to 

U plate “ tlQCOSWf 


where Uq is the maximum plate speed and u; is the frequency. The exact solution for 
this problem is given by 


u(y, 0 = «oexp 



( 6 . 1 ) 


where u is the kinematic viscosity. This is a 1-D problem in that there are no gradients 
in the direction parallel to the plate and so only rj gradient terms were used. The mesh 
extended to a height where the exact solution for the velocity is less than 5xl0 -4 «o 
and consisted of 50 points in the y-direction stretched such that the spacing at the 
plate was 0.002 of the total height. The velocity tto was set to unity, the frequency 
was set to 27r, and the square root term was set to unity. Thus 
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The physical time step At was set to 0.01. Since the vertical velocity component v is 
always zero, and there are no gradients in the x-direction, the divergence of velocity is 
always zero, and thus only one iteration is required to advance the equations one step 
in time. In fact, after one iteration, the maximum value of the right-hand side of Eq. 
(2.15) was reduced to 10 -10 . 



Fig. 8. Velocity profiles for flow over an oscillating plate. Solid line: computed 
solution; Symbols: analytical solution; o: t =8.25; A: t =8.5; +: t=8.75; x: t=9. 

The problem was run over nine full periods after which the transient solution 
had died off and the solution was fully periodic in time. In Fig. 8 velocity profiles at 
different times during the cycle are plotted for both the computed solution and the 
exact analytical solution. The computed solution is designated by the solid lines and 
the exact solution by the symbols. Excellent comparison is seen between the two. 
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6.6 One-Dimensional Channel 


The flow through an inviscid 1-D channel with an oscillating back pressure was 
calculated. This problem was also used 'by Merkle and Athavale 5 to verify their 
method. The present work used the 2-D code by setting all gradients in the y-direction 
to zero and using an equally spaced grid of 21 points in the x-direction. At the inflow 
constant total pressure was specified, and the static pressure at the outflow was set to 
be 

Pexit = PO "H Pe& 


where p e is the amplitude of the pressure oscillations and u> is the frequency. An 
exact analytical solution exists for this as long as the magnitude of the back-pressure 
oscillation p e remains small compared to the mean back pressure Pq. For this problem, 
the ratio of p e to po is taken to be 0.1. The flow was calculated for a channel length 
of unity and a mean velocity of unity. In this case the solution is given by 

i(<) = 1 — - — - (sinu>< — wcoswf — we -4 ) 


u[ 


p(x, t) — p 0 + p e sinu;< + (x — o (cos wt + wsin ut + e t ) 


( 6 . 2 ) 


1 + U> 2 

Both the pressure and velocity have exponentially decaying terms corresponding to 
the initial transient, as well as the periodic sine and cosine terms. The velocity is a 
function of time only, which is a direct effect of the incompressible continuity condition 
in a constant area duct. 


The initial conditions for the problem were taken to be the velocity and pressure 
given by Eqs. (6.2) evaluated at t=0. The problem was calculated for six different 
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frequencies ranging from 0.01 to 100.0. The physical time step was taken so that 30 
steps would be used for one period of the back-pressure oscillation. Thus 



The artificial compressibility constant /? was set to 10. Each sub-iteration in pseudo- 
time was converged until the maximum divergence of velocity at any point was less 
than 10 -6 and the maximum element of the right-hand side of Eq. (2.15) was less 
than 10 -6 . This required five to eight subiterations to obtain this convergence for all 
cases calculated. 



Fig. 9. Velocity for the flow through an inviscid channel with oscillating back pressure. 
Solid line: computed solution; X : analytical solution. 


The velocity is plotted as a function of time in Fig. 9. Shown are the computed 
values (solid line), and the exact values (x symbol) for a frequency of u> = 10. The 
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transient is seen to die out after about five periods, after which simple harmonic 
motion occurs. The velocity results compare very well with each other. For example, 
the computed mean velocity for the periodic motion is equal to 1.00035, and the exact 
mean velocity is 1. The values of the computed pressure are seen to be quite accurate 
as well. The maximum error of the computed pressure was measured for all cases ran 
and was found to average less than 10 4 , which corresponds to an error of about 0.01 
percent. This error is proably due to truncation. 

6.7 Vortex Shedding Behind a Circular Cylinder 

The flow over a circular cylinder was computed in an effort to induce vortex shed- 
ding in the wake. Computations were carried out at two different Reynolds numbers, 
105 and 200, which are discussed as separate cases. 

6.7.1 Reynolds Number 200 

The Reynolds number 200 case was run first and quantitative measurements 
of the flow field were taken. This case was run using a 60 by 100 O-grid with points 
clustered near the body such that the spacing next to the surface was 0.0035 diameters. 
The outer boundary extended to 10 diameters from the cylinder, a physical time step 
At of 0.025 was used, and the artificial compressibility constant (3 was set to 2500. 
Both the fifth- and third-order accurate convective flux differences were used. 

At each time step the subiterations were carried out until both the maximum 
divergence of velocity at any point and the maximum residual at any point were less 
than 10“ 4 . During the initial transient, 20 to 30 subiterations were required for each 


55 



time step, but this quickly decreased to an average of 5 subiterations per time step. 
Once an asymmetric wake developed, 12 to 14 subiterations were required. 

The flow was started impulsively from free-stream conditions and run until a 
periodic shedding of vortices occurred. The vortex shedding developed without any 
type of artificial disturbance. This is most probably due to the roundoff error adding 
some very small asymmetric disturbances which will trigger the asymmetric motion. 
However, the fifth-order differencing calculations started to develop an asymmetric 
wake within a nondimensional time of 50 and was completely periodic by a time of 
100, whereas the additional dissipation of the third-order scheme delayed the formation 
of a completely periodic solution until a nondimensional time of 180. 



TIME 

Fig. 10. Lift and drag coefficients versus time for flow over a circular cylinder at a 
Reynolds number of 200. 
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The lift and drag coefficients for the fifth-order calculations are plotted versus 
time in Fig. 10. The time history for the third-order calculations is quite similar to 
this plot. However the quantitative results do vary. This is shown in Table 7 where 
the values of the lift and drag coefficients for the periodic state are listed, as well 
as the Strouhal number. Also listed in this table are the results from several other 
calculations by Rosenfeld [42], Lecointe and Piquet [43], Martinez [44], Lin et al [45], 
and Thoman and Szewczyk [46]. Experimental results of Wille [47], Kovasznay [48], 
and Roshko [49] are included as well. The present fifth-order results are in general 
agreement with most of the other results. Values for the Strouhal number near 0.19 
appear to be consistent, and so the Strouhal number of the present result for the third- 
order differencing is to be too low. The difference between the third- and fifth-order 
schemes is attributed to the difference in the amount of numerical dissipation which 
is added by each scheme. 


Table 7 Lift and Drag Coefficients and Strouhal Numbers for 
Circular Cylinder Flow at Reynolds Number 200 



C D 

Cl 

St 

Present 3rd order 

1.29 ± 0.05 

± 0.75 

0.16 

Present 5th order 

1.23 ± 0.05 

± 0.65 

0.185 

Rosenfeld [42] 

1.46 ± 0.05 

± 0.69 

0.211 

Lecointe & Piquet [43] 




2nd order 

1.46 ± 0.04 

± 0.70 

0.227 

4th order 

1.58 ± 0.0035 

± 0.50 

0.194 

Martinez [44] 

1.27 ± 0.0035 



Lin et al. [45] 

1.17 



Thoman & Szewczyk [46] 

1.17 ± 0.005 



Wille [47] (exp) 

1.3 



Kovasznay [48] (exp) 



0.19 

Roshko [49] (exp) 



0.19 
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Fig. 11. Streamlines for flow over a circular cylinder at a Reynolds number of 200 at 
various times during the vortex shedding cycle for the third-order differencing scheme. 

In Fig. 11, the streamlines at various stages during one period are plotted for the 
third-order calculations. The first plot shows the flow when the drag is at a minimum 
and the lift is zero. It can be seen that a new vortex is forming on the top half, while 
the low-pressure center of the previous vortex has pulled away from the body causing 
the drop in drag. The second plot shows that the top vortex has now extended itself 
across the entire back side of the body, causing this point in time to have a maximum 
in drag and lift. The next two plots correspond to another minimum in drag with 
zero lift, and a maximum in drag with a minimum in lift, respectively. These two are 
mirror images of the first two plots. Finally, the last plot is nearly identical to the first 
one, and is only different because it is taken at a time of 6 after the first plot, whereas 
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the actual period for this flow is 6.25. The main qualitative flow features shown here 
are the same as are seen in the fifth-order calculations. 


6.7.2 Reynolds Number 105 

The code was next run for a circular cylinder at a Reynolds number of 105 in an 
effort to see how well the Karman vortex street is captured by the upwind differencing 
scheme. The grid dimensions were increased to 120 by 240 points in the radial and 
circumferential directions, respectively, and the outer boundary was extended to 25 
diameters away from the cylinder. A physical time step of 0.05 was used, and the value 
of / 3 was set to 1100. The fifth-order convective flux differences were used. In order to 
reduce the computing time required to see the vortex shedding, an artificial velocity 
disturbance along one side of the cylinder was added. The surface velocity on one 
side was prescribed to be sinusoidal in time with a non-dimensional frequency of 0.16 
and a magnitude of one quarter of the freestream velocity. This triggered shedding 
right away, and the disturbance was turned off after a non-dimensional time of 15.7, 
after which time the vortex shedding continued and a completely periodic solution 
was obtained within several cycles. With this finer grid, an average of 30 subiterations 
were required at each time step. The computing time requirements for this case were 
about 45 minutes per periodic cycle. 

Streaklines were computed for this periodic solution at a non-dimensional time 
of 58 and are shown in Fig. 12. The lower half of this figure shows an experimental 
picture of the same flow conditions by Sadatoshi Taneda reproduced from Van Dyke 
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Fig. 12. Comparison between computational and experimental streaklines for the flow 
over a circular cylinder at a Reynolds number of 105. 

[50] . The streaklines in the experiment are shown by electrolytic precipitation in water. 
The vortex structures are seen to be very similar. 

6.8 Square Duct With 00° Bend 

The flow through a square duct with a 90° bend was used as a steady-state 3-D 
test case. This particular geometry was studied experimentally by Humphrey et al 

[51] which enables comparisons to be made with the current computed results. Four 
different grids were used whose dimensions are 31 x 11 x 11, 41 x 21 x 21, 51 x 31 x 31, 
and 61 x 41 x 41. The problem was non-dimensionalized using the side of the square 
cross-section as the unit length, and the average inflow velocity as the unit velocity. 
The Reynolds number was 790 based on the unit length and velocity, and the artificial 
compressibility parameter /3 was set to one for all cases. The grid for the 31 x 11 x 
11 case is shown in Fig. 13. The straight inflow section before the bend was set to a 
length of five and the outflow section downstream of the bend was also set to a length 
of five. The radius of curvature of the inner wall in the curved section was 1.8 units 
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Fig. 13. Geometry and grid (31 x 11 x 11) for computation of flow through a square 
duct with a 90° bend. 

in length. The inflow velocity profile was prescribed. to be that of a fully developed 
laminar straight square duct as given by White [52]. 

The convergence history for this steady-state problem is shown in Fig. 14. The 
log of the maximum residual over all the grid points is plotted versus iteration number 
for all four grid cases. The convergence is shown to be very fast. Although not shown 
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Iteration Number 

Fig. 14. Maximum residual over all points versus iteration number showing conver- 
gence history for the flow through a square duct with a 90° bend. 


in the figure, machine zero is reached in less than 300 iterations for all four cases. The 
solution is considered converged if the maximum residual has converged at least four 
orders of magnitude and the maximum divergence of velocity is less then 10“ 4 . This 
is obtained in less than 115 iterations for all cases. The computing times required for 
convergence for the four cases are 32, 196, 550, and 1067 sec, respectively. 
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The computed results are compared to the experimental results of Humphrey et 
al. [51] in Fig. 15. Shown are the longitudinal velocity profiles at various streamwise 
stations for two different cross-flow locations. In Fig. 15a, the velocity profiles are 
taken at z=0.25, that is half way between the x-y plane wall and the x-y symmetry 
plane. Th- second location, shown in Fig. 15b, is from the x-y plane at z=0.5, or the 
x-y symmetry plane. In each of these figures, the profiles are shown at x=0 (the inflow 
boundary), at x=2.5 (half-way between the inflow boundary and the start of the curved 
section), and at four positions in the curved section corresponding to 6 equal to 0, 30, 
60, and 90.° The symbols represent the experimental results and the lines represent 
the computed solutions. The computations for the two finest grid cases are seen to be 
in very good agreement, indicating that these represent a grid-independent solution. 
Good comparison is seen between the computation and the experiment, particularly 
for the first four streamwise stations. However, at the latter part of the bend, the 
fine-grid computations begin to predict two relative minima in the velocity profiles, 
which is not evident in the experimental results. It is seen that the computations are 
able to predict the maximum longitudinal velocity very well. 


Figure 16 shows the cross-stream velocities for the 51 x 31 x 31 grid case at the 
streamwise points of 6 = 30°, 9 = 60°, and 6 = 90°. This figures shows how a pair of 
secondary vortices are generated by the large values of static pressure on the outside 
wall, which arises when the flow starts to negotiate the sharp bend. The center of these 
vortices is seen to move towards the inside wall between the 6 = 30° station and the 
$ = 60° position. The vortices tend to center again further downstream, and at the 
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R R 

Fig. 15. Streamwise-velocity profiles at the strearawise stations given by x=0, x=2.5, 
0 = O°,0 = 30 °, 0 = 60°, and 9 = 90°. Solid line = 31 x 11 x 11 grid, dash line = 41 x 
21 x 21 grid, dotted line = 51 x 31 x 31 grid, chain-dot line = 61 x 41 x 41 grid, and 
o = experiment [51]. 


same time a secondary pair of vortices are seen to appear. This agrees qualitatively 
with the observations of Humphrey et al. [51]. 

Figure 17 shows velocity vectors for the 51 x 31 x 31 grid case at the k=2 grid 
plane next to the outer-radius wall between the streamwise positions of 9 = 0 and 9 = 
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r = 2.8 


r = 1.8 1! 



9 = 60 ° 



Fig. 16. Cross-sectional velocities for the curved square duct at the streamwise stations 
given by 6 = 30 °, 0 = 60°, and 9 = 90° for 51 x 31 x 31 grid. 


60.° These show the presence of a recirculation region in the corners of the duct 
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Fig. 17. Velocity vectors next to outer-radius wall showing reverse flow in curved 
square duct for 51 x 31 x 31 grid. 


appearing between 6 = 0 and 6 = 35.° In their experiment, Humphrey et al. [51] 
reported recirculating flow in the corners between approximately 6 = 0 and 6 = 25.° 


6.0 Artificial Heart Flow 


The present flow solver has been used to compute the flow inside an artificial 
heart. The artificial heart was designed by Penn State University and is being studied 
experimentally by Tarbell et al. [53]. The purpose of the current calculations is to 
demonstrate and analyze the present capability to compute a time-accurate incom- 
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pressible flow through a complex internal device with moving-boundaries. The initial 
calculations in this effort are presented here. Since a number of primarily geometrical 
differences exist between these initial calculations and the actual artificial heart exper- 
iment, which are detailed below, only the simplest of qualitative comparison between 
computation and experiment can be made. Further work will attempt to remove these 
differences so that the actual artificial heart geometry can be accurately modeled. 



Outflow Valve Inflow Valve 

Fig. 18. Artificial-heart geometry showing valve openings. 


The geometry used for the current model is depicted in Fig. 18. The heart 
is composed of a cylindrical chamber with two openings on the side for valves. The 
pumping action is provided by a piston surface which moves up and down inside the 
chamber. The diameter of the piston is 7.4 cm, with a stroke length of 2.54 cm. The 
problem was nondimensionalized with a unit length of 2.54 cm and a unit velocity of 
40 cm/sec. The actual artificial heart has cylindrical tubes extending out of each of the 
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side valve openings. These contain tilting flat disks which open and close to act as the 
valves. In the computational model these valves are not modeled, instead the boundary 
conditions at the side openings are specified to instantaneously open and close at the 
right moment. This simplification allows a single zone to be used to model the flow 
inside the chamber. With the addition of more zones, however, it will be possible to 
model the tilting disk valves as well. Additional simplifications made in the present 
computational model include the movement of the piston. In the actual device the 
piston moves through the entire chamber volume, including across the majority of the 
valve openings. Because of grid-generation problems, moving the piston past the valve 
openings would become quite difficult, and so the piston is restricted to move only 
in the volume beneath the valve openings. The bottom most position of the piston 
is extended lower than normal so that the entire piston displacement volume is close 
to that in the actual artificial heart. The flow is assumed laminar, and the Reynolds 
number based on the the unit length and velocity is set to 100. In the actual heart the 
Reynolds number is about 600, and regions of the flow are turbulent. Finally, the fluid 
is assumed to be Newtonian. This corresponds to the experiment of Tarbell et al. [63] 
who used a water and glycerin fluid whose viscosity is nearly the same as blood, about 
3.5 centipoise, but unlike blood can be simulated by a Newtonian fluid assumption. 

Inside the heart an H-H grid topology is used with dimensions of 39 x 39 x 51 
Figure 19 shows the grid on the unwrapped surface of the side of the heart chamber. 
Grid lines were placed along the lines of the valve openings to make boundary condition 
implementation straight forward. This surface grid was generated by dividing the 
surface into several zones and using a biharmonic grid generator [54] in each zone. 
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Fig. 19. Unwrapped grid used along side of the artificial heart near the region of the 
valve openings. 

The same grid generator was also used to generate an H grid on the piston surface 
and on the top surface. To fill in the interior points, an algebraic solver coupled with 
an elliptic smoother was used. As the piston moved up and down inside the chamber, 
the grid points below the valve openings were compressed and expanded, respectively. 
Thus a new grid was generated at each time step. 

The flow was computed using a time step At of 0.025, and a /? of 500. The 
piston moved with a constant nondimensionalized velocity of ±0.2 between its top and 
bottom positions, thus requiring 200 physical time steps for one period of the piston’s 
motion. During each time step, the subiterations were carried out until the maximum 
residual dropped below 10 -3 or until a maximum of 20 subiterations were used. During 
most of the piston’s cycle only 12-15 subiterations were required, but when the piston 
was changing directions, it did not completely converge in 20 subiterations. This did 
not cause any stability problems, yet it remains to be seen what effect this might have 
had on the time accuracy. The computing time required for each period of the piston’s 
motion was approximately four hrs. The computations were rim for four periods during 
which particle paths were computed after being released near the inflow valve. Figure 



20b. Experiment 


Fig. 20. Incoming particle traces from computations and picture of experimental 
results as the piston nears the bottom position. 


20a shows some of these particle traces as the piston nears its bottom position. Two 
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distinct vortices are seen to have formed from the flow separating as it enters through 
the inflow valve. In Fig. 20b, an experimental photograph (Tarbell, J. M.: private 
communication, 1988) shows bubbles entering the inflow valve as the piston nears its 
bottom position. A similar two-vortex system is seen to form here. Figures 21 and 22 
show velocity vectors during the inflow phase in planes passing through the center of 
the inflow valve. The first of these shows a top view of vectors in a plane parallel to 
the piston, while Fig. 22 shows a side view of vectors in a plane perpendicular to the 
piston. These figures portray how complicated the vortical structure of this flow is. 
Figure 21 again shows the presence of two vortices formed as the incoming flow forms 
a jet. Figure 22 shows also how the flow separates underneath the valve opening, but 
what it does not show is that the flow there is strongly 3-D with the velocity vectors 
next to the left wall underneath the valve pointing into the paper. Also in the figure 
is seen the presence of additional vortices against the back wall opposite the valve 
opening. 
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Fig. 21. Top view of velocity vectors in plane through center of inflow valve showing 
incoming fluid. 
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Fig. 22. Side view of velocity vectors in plane through center of inflow valve showing 
incoming fluid. 
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CHAPTER SEVEN 


CONCLUSIONS 
7.1 Summary 

An algorithm for computing both steady-state and time-accurate solutions to the 
incompressible Navier-Stokes equations has been presented. This approach has been 
shown to be capable of computing flow about 3-D realistic geometries with moving 
boundaries. The method of artificial compressibility allows the equations to be solved 
as a hyperbolic system in pseudo-time. This requires the solution of a steady-state 
problem at each physical time step for the time-accurate formulation. The use of up- 
wind differencing makes the system of numerical equations more diagonally dominant 
than a central-differencing scheme, leading to significantly faster convergence and far 
less computing time requirements. This is possible partially because of the use of a 
nonfactored implicit Gauss-Seidal line-relaxation scheme, making it possible to run 
the upwind-differencing scheme at very large time steps. Additional advantages of 
the current code include the smaller number of numerical parameters which must be 
specified. For this code, only the artificial compressibility (3, and the number of line- 
relaxation sweeps need to be specified. These parameters are generally independent, 
that is, changing the value of one will not dramatically change the optimum value of 
the other. In contrast, an artificial compressibility code which uses central differencing 
and approximate factorization requires the specification of three dependent parame- 
ters: the artificial compressibility constant /3, the artificial dissipation parameters, and 
the time step size. 
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The results showed good comparison with experimental results and with ana- 
lytical solutions for the eight laminar flow test cases presented. Computing the flow 
through an artificial heart shows the capability of the code to simulate complicated 
internal flows with moving boundaries within a reasonable amount of computing time. 

The choice of the artificial compressibility parameter ft has been found to be 
very important. Its value can affect the convergence rate dramatically. The optimum 
values of ft were found in the present work using numerical convergence tests. Two 
general dependencies were observed in this process: that the optimum ft is inversely 
proportional to the Reynolds number of the flow problem; and that ft is proportional 
to the grid density. 

7.2 Future Work 

Further advances in the convergence speed of the algorithm will still be very use- 
ful in increasing the usefulness of this code as a design tool. Several approaches could 
lead to some improvement in this area. The first and perhaps easiest to implement 
would be the use of a multigrid acceleration scheme. Some less obvious improve- 
ments could be sought after by studying the possibility of diagonalizing the left-hand 
side tridiagonal-system of equations which are solved for each separate line in the 
line-relaxation process. This could reduce the work of solving 4x4 blocks (in 3-D) to 
solving 4 independent scalar equations. 

Another concern is the large memory resources required by the current implemen- 
tation of the line-relaxation scheme in which all of the left-hand side terms are formed 
at once and stored. This implementation is made easier because the incompressible 
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Navier-Stokes have four dependent variables in 3-D, in contrast to the compressible 
equations which have five. This large memory usage makes this code useful for 3-D 
calculations on only a few of the supercomputers in use today. Some work will have 
to be done to test the amount of additional time required to compute the implicit-side 
terms as they are required instead of all at once. Testing of the artificial compressibil- 
ity method using upwind differencing with several other types of implicit solvers would 
be very useful in further determining this algorithm’s usefulness on supercomputers 
where memory is not as abundant. Some of these schemes could include approximate 
factorization and an LU factored scheme, such as the one by Yoon and Jameson [55]. 
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APPENDIX A 


VISCOUS FLUXES 


This appendix presents the differential formulas for the viscous fluxes in general- 
ized coordinates in both 2-D and 3-D for the following sets of conditions: nonconstant 
viscosity on a nonorthogonal mesh; nonconstant viscosity on an orthogonal mesh; con- 
stant viscosity on a nonorthogonal mesh; and constant viscosity on an orthogonal 
mesh. In the following equations, the velocity gradients in the viscous fluxes were 

written as 

du 

— = u ( , etc. 

and the metrics of the transformation have been represented with 


Three-Dimensional Formulation 


Nonconstant Viscosity, Nonorthogonal Mesh 
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Nonconstant Viscosity, Orthogonal Mesh 
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Constant Viscosity, Nonorthogonal Mesh 
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Two-Dimensional Formulation 


Nonconstant Viscosity, Nonorthogonal Mesh 
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Nonconst ant Viscosity, Orthogonal Mesh 
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APPENDIX B 


EIGENSYSTEM OF THE CONVECTIVE FLUX JACOBIAN 

3-D System of Equations 

To form the delta fluxes used in the flux- differ erce splitting scheme, the eigensys- 
tem of the convective flux Jacobian is needed. For the equations in 3-D a generalized 
flux vector is given by 


Ei = 
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L k t w + k z p + w QJ 


(B. 1 ) 


where E* = E,F,G for i = 1,2, and 3 respectively, and the normalized metrics are 
represented with 
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The Jacobian matrix Ai = ^ of the flux vector in Eq. (B.l) is given by 

's x p kyp k z p 

7 I A-z ^2“ T Q + kt kyU k z 1l /n M 

Ai — 7 l „ u j n _l u u 


• 0 

k x P 

kyp 

k z p 

K 

k x ii 4 " Q kf 

kyU 

k z u 

ky 

k x v 

kyV + Q + kt 

k z v 

K 

k x w 

k y w 

k z w + Q + fc, 


82 


A similarity transform for the Jacobian matrix is introduced 


Ai=XiAiX7' 


where 

Ai = diag[Ai, A 2 , A 3 , A 4 ] 

Ai = Q + kf 

A2 = Q + kt (i?. 3 ) 
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where 
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2-D System of Equations 


A generalized flux vector is given by 


m 

Ei = I k x p + uQ + k t u ( B.7 ) 
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where Ei = E, F for i = 1,2 respectively, and the normalized metrics are represented 
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A similarity transform for the Jacobian matrix is introduced 


Ai = XiAiX- 


where 

A i = diag[Ai, A 2 ,A 3 ] 

Aj = Q + kt 

1, 

A2 — Q + c + -k t 

A 3 = Q - C + ^ kt 

and where c is the scaled artificial speed of sound given by 
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c = yW + + tj) 


(5.10) 


The matrix of the right eigenvectors is given by 


Xi = 


1 

2 M^W) 


0 /3(c 2 - \k\) 

-2j3ck y (u\ 2 + 0k x )(c + ifcf) 
2/3ck x (vA 2 + (3k y )(c + |fc t ) 


~/?(c 2 - jA: 2 ) 

(ixA 3 + /?A; x )(c- 
(uA 3 +/?fcj,)(c- | A:*) 

(5.11) 


and its inverse is given by 



k y u — k x v —v\ 1 —/3ky u\i + (3k x 
— A 3 (3k x (3k y 

— A 2 /?A: x . /5fcj, 


(5.12) 
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